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ABSTRACT 

We present results based on an implemen tation of the Go dunov Smoothed Particle Hydrody- 
namics (GSPH), originally developed by llnutsukal (l2002h . in the GADGET-3 hydrodynamic 
code. We first review the derivation of the GSPH discretization of the equations of moment 
and energy conservation, starting from the convolution of these equations with the interpolat- 
ing kernel. The two most important aspects of the numerical implementation of these equa- 
tions are (a) the appearance of fluid velocity and pressure obtained from the solution of the 
Riemann problem between each pair of particles, and (b) the absence of an artificial viscosity 
term. We carry out three different controlled hydrodynamic al three-dimensional tests, namely 
the Sod shock tube, the development of Kelvin-Helmholtz instabilities in a shear flow test, 
and the "blob" test describing the evolution of a cold cloud moving against a hot wind. 

The resu l ts of o ur tests confirm and extend in a number of aspects those recently obtained 
by ICha et all (|2010|) : (i) GSPH provides a much improved description of contact disconti- 
nuities, with respect to SPH, thus avoiding the appearance of spurious pressure forces; (ii) 
GSPH is able to follow the development of gas-dynamical instabilities, such as the Kevin- 
Hehnholtz and the Rayleigh-Taylor ones; (iii) as a result, GSPH describes the development of 
curl structures in the shear-flow test and the dissolution of the cold cloud in the "blob" test. 

Besides comparing the results of GSPH with those from standard SPH implementations, 
we also discuss in detail the effect on the performances of GSPH of changing different as- 
pects of its implementation: choice of the number of neighbours, accuracy of the interpolation 
procedure to locate the interface between two fluid elements (particles) for the solution of 
the Riemann problem, order of the reconstruction for the assignment of variables at the in- 
terface, choice of the limiter to prevent oscillations of interpolated quantities in the solution 
of the Riemann Problem. The results of our tests demonstrate that GSPH is in fact a highly 
promising hydrodynamic scheme, also to be coupled to an N-body solver, for astrophysical 
and cosmological applications. 

Key words: Hydrodynamics - instabilities - turbulence - methods: numerical - galaxies: 
formation. 



1 INTRODUCTION 

Lagrangian hydrodynamic methods find a natural field of applica- 
tion in astrophysics and cosmology. Thanks to their intrinsically 
adaptive nature, they are well suited to deal with large ranges 
of scales and densities, as well as with the complex geometries 
of the typical problems involved in astrophysics. Among such 
Lagrangian methods. Smoothed Particl e Hydrodynamics (SPH; 
iGingold & Monaghan|[l977l ; lLucvlll977h has b een, and currently 
is, by far the most widely used scheme (e.g., iMonaghMl l2005l : 
iRosswodlioO^ ; ISpringelll2010bl , for recent reviews on SPH and on 



its applications to astrophysical problems). The fluid representa- 
tion by particles that move with the flow makes SPH providing a 
solution of the Euler equation for an inviscid fluid in Lagrangian 
coordinates. Remarkable advantages of this representation are its 
intrinsic Galilean invariance and the possibility to easily couple it 
to N-body solvers describing the dynamics of self-gravitating flu- 
ids, or to a fluid feeling an external gravitational potential. 

Despite such advantages, a number of intrinsic limitations of 
SPH have been recognised. Historically, the first one is related to 
the difficulty that SPH has in capturing shocks and contact disconti- 
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nuities. At shock fronts, the Rankine-Hugoniot jumping conditions 
predict specific entropy of the fluid to increase through the con- 
version of mechanical energy into internal energy, thus implying 
that an inviscid description of the fluid can not be valid any longer. 
To deal with this problem, one needs to introduce in SPH an arti- 
ficial viscosity to dis sipate local velocity differen c es and convert 
them into heat (e.g., Monaghan & Gingold [19831 : lBalsaralll993 : 
lMonaghanlll997l : lLombardi et al.l ll999^). The typical viscosity re- 
quired for this is larger than the natural viscosity of the fluid. This 
leads to the danger that it may also provide spurious effects away 
from shock regions, e.g. causing an unphysical damping of turbu- 
lent motions or spurious angular momentum transport in differen- 
tially rotating discs. Therefore, any scheme of artificial viscosity 
need to be tuned so as to be localised as much as possible at shock 
regions. Indeed, attempts have been devoted to design schemes 
in which artificial viscosity is reduced or eve n eliminated away 
from shocks (e.g. , Morris & M onaghan 1997: IPolag et al.ll2005l : 
ICullen & Dehne j|2010l) . 

Another, possibly more serious limitation, of SPH lies in its 
difficulty to correctly describe fluid instabilities and mixing at 
the boundaries between different fluid phases. The main reasons 
for this are due to the limitations of SPH in computing gradients 
across discontinuities and to the intrinsic lack of diffusivity of SPH, 
which ma kes entropy to be conserved within the kernel scale (see 
iRead et a l. 2010, and references therein, for a detailed discussion ). 
This limitation of SPH has been highlighted by Agertz et al.' llQQl). 
These authors caiTied out a comparison between different SPH and 
Eulerian grid codes, with the aim of comparing their relative capa- 
bility of describing the onset of fluidodynamical instabilities in spe- 
cific test cases. One of the main results of this analysis was that the 
incorrect description of contact discontinuities provided by SPH 
causes a sort of spurious surface tension to appear, which prevents 
the development of Kelvin-Helmholtz (KH) and Rayleigh-Taylor 
(RT) instabilities. 

Stimulated by this analysis, a number of aut hors propose d dif- 
ferent approaches to improve SPH performances. IPricafcOOSh sug- 
gested that discontinuity of entropy at fluid interfaces, combined 
with a smooth variation of density, causes the appearance of "pres- 
sure blips" as a consequence of a spurious surface tension. To solve 
this problem, he introduces a thermal conduction term which acts 
as a diffusive term in the energy equation, thus improving the capa- 
bility of SPH of describing mixing and the development of instabil- 
ities (see also|M erlin et al. 2010). A similar argument was also pre- 
sented by Wadsley et al. ( 2008), who resorted t o a sub-grid model 
of turbulence as the source of thermal diffusion. iRead et al] 1 I2OIOI) 
adopted a different approach. They employed a different kernel, a 
much larger number of neighbours, and a modified density estima- 
tion formula to control the errors in the estimates of gradients, so 
as to obtain a better representation of mixing of different phases in 
shear layers. 

Wel l befor e this vivid debate on the limitations of SPH started, 
llnutsukal j2002l) (102 hereafter) proposed a novel approach to La- 
grangian hydrodynamics. This approach was based on two main 
considerations. First, the standard SPH approach is based on as- 
suming that coarse-grained thermodynamic quantities (i.e. density, 
pressure, entropy), assigned at the particle positions by kernel con- 
volution, can be evolved through the equation of fluido-dynamics. 
However, strictly speaking, equation of fluido-dynamics describe 
the evolution of micro-physical (i.e. non coarse-grained) quan- 
tities. Therefore, a self-consistent particle description of fluido- 
dynamics would require the coarse-graining procedure to be ap- 
plied to the equations evolving micro-physical variables, rather 



than to the micro-physical variables themselves. The two opera- 
tions, i.e. coarse-graining and dynamical evolution, in general do 
not commute. It is the violation of this commutation that originates 
the approximation in the fluido-dynamical description provided by 
SPH. Second, in deriving the implementation, through particle de- 
scription, of the equations of evolution, 102 devoted special atten- 
tion to keep the order of spatial accuracy (i.e., in kernel smooth- 
ing length) fixed to second-order. The natural way of implement- 
ing this SPH formulation is by computing the exchange of hydro- 
dynamic forces and momenta between pairs of particles by using 
a Riemann sol ver, analogous t o the grid-based second-order Go- 
dunov scheme jvan Leeiill979h . Therefore, the resulting Godunov 
Smoothed Particle Hydrodynamics (GSPH hereafter) method has 
the extra benefit of not requiring the introduction of an artificial 
viscosity term, since shocks are now naturally described as so- 
lutions of the Riemann problem (RP hereafter). 102 showed with 
one-dimensional tests that the mixing associated to the energy and 
momentum exchange naturally prevents t he formation of "pressur e 
blips" at conta ct discontinuities (see also lCha & Whitworthll2003h . 
More recentlv. lcha et aklfcoiOl) showed that GSPH also provides a 
much improved description of the development of KH instabilities 
through two-dimensional tests. 

In this paper we present results based on o ur implementa - 
tion of GSPH on the GADGET-3 simulation code ISpringelll2005h . 
Using standard three-dimensional hydrodynamic tests, we aim at 
demonstrating the capability of GSPH of describing mixing and 
development of fluidodynamical instabilities at interfaces. In our 
analysis we will pay special attention to i) highlight the fundamen- 
tal differences with respect to the standard SPH approach; ii) dis- 
cuss the effect of changing relevant aspects of the implementation 
of the Riemann solver^ 

We note that Springell ( l2010al) recently proposed a novel 
particle-based method, in which particle positions are used to con- 
struct an unstructured moving mesh, according to the Voronoi tes- 
sellation. Similarly to the GSPH, also in this case fluid equations 
are solved with the finite volume Eulerian Godunov scheme, with 
fluxes across the boundaries of the Voronoi poliedra provided by 
the solution of the Riemann problem between particle pairs. While 
this scheme and GSPH appear to be similar in spirit, there are a 
number of fundamental differences. The most important is proba- 
bly represented by the fact that, while GSPH always refers to the 
volumes as defined by the k ernel smoothing length, the scheme 
proposed bv lSpringeil ( I2010al) uses the partitioning of the computa- 
tional domain given by the Voronoi tessellation. 

The scheme of our paper is as follows. Section 2 gives the ba- 
sis of our GSPH formulation. In Section 3 we describe the details 
of our implementation of GSPH in GADGET. In Section 4 we de- 
scribe our results for the hydrodynamical test of the shock tube, the 
shear flows and of the "blob" test. In Section 5 we draw our con- 
clusions. The Appendix contains the expression for interpolating 
volumes and for the position of the interface. 



2 BASICS OF THE GODUNOV SPH 

In the following we provide a short description of the approach at 
the basis of t he Godunov SP H (GSPH) method, while we refer to 
the paper bv llnutsukal ( l2002h for a more complete formal deriva- 
tion. 

Let us introduce the convolution of a physical function /(x) 
with the kernel function. 
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/(x')VF(x-x';;i)dx', 



(1) 



where, as usual, h denotes the kernel size, that for the moment we 
assume to be spatially constant, and the integral is performed over 
the whole spatial domain. Starting from the above expression, it is 
easy to demonstrate that (V/) = V(/). 

Defining the density field at the position x as given by the 
sum of the kernel contributions at particles positions Xj, p(x) — 
. rrij W{x — ■Kj-, h), we have then the identities 



E 



^W{^^^';h) ; = Vm,V 



W{x-xf,h) 



.(2) 



Using the first of the two above identities, the expression of the 
kernel convolution can be cast in the form 

/, = (/)(xO = / Vm, 4^1^(x'-x,;/i)Vy(x'-x,;/i)dx'.(3) 

We note that the expression for the kernel convolution in the stan- 
dard SPH approach is recovered under the approximation VK(x — 
x';/i) =5b(x-x'). 

To derive the equation of evolution for particles, we start from 
the kernel convolution of the equation of motion, 

I ^XM (x-x'; h)d^ = ^ [ VK(x-x'; h) dx , (4) 

J dt J p{x) 

where v(x) is the velocity field and P(x) the pressure field of the 
fluid. Since Xi = J W(x — x'; /i) dx describes the motion of 
the i-th particle position, integrating by part the r.h.s. of the above 
equation and using the first identity in Eq.(|2j, we obtain the follow- 
ing expression for the equation of motion: 



m-iXi 



P(x) 

p(x)2 



[^^ ~ dj] WrWj dx . 



(5) 



where we introduced the notations di = d/d'x.i and Wi = W^(x — 
Xi; h). 

As for the energy equation, the coarse-grained representation 
for the evolution of the specific internal energy m(x) can be written 

as 



dit(x) 



dt 
P(x) 
P(x) 



VK(x-x';ft)dx : 



[V-v] W(x-x';/!.)dx. 



(6) 



After rearranging the r.h.s. of the above equation and using the ap- 
proximation 



V ■ VP(x) 

P(x) 
Xi • VP(x 

P(x) 



W^(ix-x'|;ft)dx = 

■ W(|x-x'|;/0dx + O(;i^), 



(7) 



the equation for the evolution of the internal energy of the i-th par- 
ticle can be written as 



{v - Xi) ■ {d, - d,)WiWj 



(8) 



We point out that Eqs.l|5j and ^ replace in the GPSH ap- 
proach the corresponding equations of moti on and energy of the 
standard SPH approach. As demonstrated bv llnutsukal ( l2002h . they 



provide a description of the coarse-grained equations of fluido- 
dynamics of Eqs. Q and which is second-order accurate in 
spatial resolution (i.e. in kernel smoothing length). 



3 IMPLEMENTATION 

In this section we describe the numerical implementation of the 
GSPH equations of evolution. We remind that we implemented 
these equations in the massively-parallel GADGET-3 simula- 
tion code. GADGET-3 is a N-body/hydrodynamic co de in which 
entropy-conserving SPH dSpringel & Hemquisd l2002l) is coupled 
to a TreePM N-body solver to describe gravity. The code has 
fully adaptive time-stepping. Domain decomposition is carried out 
by using a space-filling Peano-Hilbert curve, which is split into 
segments assigned to different computing units. In this respect, 
GADGET-3 represents a substa ntial improveme nt with respect to 
the previous GADGET-2 version ISprineelldZOOSh , in that disjointed 
segments of the Peano-Hilbert curve can be assigned to a single 
computing unit, so as to achieve an optimal work-load balance. Be- 
sides the reference entropy conserving SPH formulation, GADGET 
also includes a switch to standard energy conserving SPH. Our im- 
plementation of GSPH replaces the equations of standard SPH. 



3.1 The GSPH equations 

102 described how to evaluate the spatial integrals appearing in 
Eqs. {5)1 and {8)1, under the assumption of Gaussian kernel. 



W{x,h) 



■d/2 -x^/h^ 

' e ' 



(9) 



where d is the number of dimensions. To this purpose, for a given 
pair of particles having coordinates Xi and Xj, one defines the s- 
axis, which is parallel to the direction of the separation vector x^ — 
Xj, with origin at (xi — Xj)/2. We also denote with Si and Sj 
the components of the Xi and Xj vectors along the s-axis, so that 
Asij — Si — Sj = \xi — Xj|. Defining then the specific volume 
occupied by a fluid element of density p(x) as V^(x) — l/p(x), 
then its gradient is Vl^(x) — — V(x — x^; /i)/p^(x). It is 

then possible to demonstrate that after expanding l/p^ to the first 
order in the direction parallel to the vector Xi — Xj , the equations 
of evolution can be rewritten in the form 



Axj 
At 



-2^mjP*K^j(/i)9,I¥(x, 



■V2h), 



(10) 



;\/2/i).(ll) 



In the above equations A indicates finite difference of each vari- 
able, Xi is the time-centred velocity of the i-th particle, while P* 
and V* are provided by the solution of the Riemann problem be- 
tween the i-th and the j-th particles. The use of P* instead of P(x) 
is justified using a linear interpolation for P(x), and evaluating it 
at the position s* ^ (See the Appendix for details). Furthermore, the 
quantity V^j accounts for the expansion of the term and can 
be expressed through the kernel convolution according to the rela- 
tion 



(12) 



/9-^(x) W{yc - XO W{X - Xj) dx : 

VK(x, -x,;\/2/i). 
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We point out that the possibility of factoring out in the above equa- 
tion the dependence on the separation vector of the particle pair, 
and the \/2 factor appearing in front of the smoothing length in the 
r.h.s. stem from the assumption of Gaussian kernel. 

We provide in Appendix the expressions for the position of 
the interface and for the Vi^j quantities in the case of linear and 
of cubic spline interpolation for the Vijx) in the s-coordinate. In 
the following we will use the more accurate cubic spline in our 
reference GSPH formulation. We will also show the effect of using 
a linear interpolation for volumes, in the shear flow and blob tests. 

We also point out that the derivation of the GSPH equations of 
evolution l llOt and i ll It have been derived by assuming a constant 
value for the kernel smoothing length h. In the case of adaptive 
smoothing, the formal derivation that leads to Eqs.(|5) and (jSj can 
be repeated by just replacing h with /i(x). However, for a spatially 
varying smoothing length the convolution integrals leading to the 
GSPH equations JlOb and (TTJ can not be done analytically. In this 
case, 102 made the ansatz that hi — h{xi) has to be used for half 
of the integration volume and hj = h{xj) for the other half. This 
leads to the following expressions for the GSPH evolution equa- 
tions in the presence of adaptive smoothing: 



At 



Vi,Jih,)^^WiK^ 



■,V2hr) 



■,V2hi 



(13) 



At 



mjP*[v* — Xi] 



V^^,{hi)diW{xi-Xj-V2h,) 



+ Vi'^,{h/)^^W{x,-XJ■V2hJ) 



(14) 



In the following we will estimate the local value of /i(x) by as- 
suming that the kernel contains a fixed number A'^„cigh of particles. 
The Gaussian kernel has not a compact support, thus implying that 
each particle should in principle interact with all the other parti- 
cles. This would result in an extremely expensive calculation, with 
a time cost tcpu oc Np^^^^, scaling with the square of the particle 
number. To avoid this, we simply truncate the Gaussian kernel at a 
distance r — 3h, the neglected contribution of the kernel being of 
the order of 10"^ 

We also discuss the effect of varying the number of neighbours 
in the case of the KH test (see Section [4~2l below). We note that this 
criterion to choose the number of neighbours is different from that 
adopted by 102, which is instead based on the requirement that the 
resulting /i(x) is not varying much within the neighbourhood of 
each particle. 

Furthermore, for each particle i the sums appearing in 
Eqs. l ll3t and il4\ must be performed over all neighbours within 
a distance \/2hi or \pihj. This further increases the number of 
neighbours thus incre asing the comp utational cost. 

As discussed by ICha et al.l ( I2OIO.) . the standard SPH evalu- 
ation of gas density can generate an unphysical repulsive force 
in some particular particle distributions, for example a non uni- 
form one. This problem is prevented by defining density as an even 
function for the exchange of particle positions. Therefore, we sym- 
metrize our density estimate with respect to the Gaussian kernel, 
according to 



p(xi; 



E 



(15) 



where = [[W{xj - x,; y^h) + W{y^^ - x,; ^/ij)] /2. 
We use this estimate of the gas density in all of our GSPH formu- 
lations, unless otherwise specified. 

We implement Eqs. (T3) and (Hi in the GADGET-3 code. Be- 
sides avoiding the need of introducing the artificial viscosity term 
in the equation of motion, another fundamental difference between 
the GSPH and the standard SPH evolution equations lies in the fact 
that velocity and pressure terms associated to each pair of particles 
are replaced by the corresponding RP solutions, evaluated at the in- 
terface position. This causes the appearance in the energy equation 
of the term between square brackets in the r.h.s., which effectively 
represents a mixing term for internal energy. This is inherently dif- 
ferent from SPH, which instead provides a strictly non-diffusive 
description of the evolution of internal energy. 

However, it is worth pointing out that the above GSPH evo- 
lution equations can not be obtained by simply replacing pressure 
and velocity terms in the standard SPH fomulation, with the values 
provided by the solution of the Riemann Problem between i-th and 
j-th particle. A further fundamental point of difference lies in the 
interpolating volumes Vi,j. These terms account for the fact that 
GSPH equations are directly derived fro m the c onvolution of the 
energy and momentum equations. Indeed lcha & Whitworth ( 20031) 
introduced a variant of these equations, which neglects the convo- 
lution integrals and, therefore, are equivalent to replacing the inter- 
polating volumes with the values oiXj (? computed at the positions 
of the i-th and j-th particle: 



Axi \ - . 

AT = E'-^^ 



+ 



9iW(xi — Xj; /ij 



(16) 



Aui 
'At 



^mj P*[v* -Xi 



diW{xi ~ •x.j;hi 



-diW{^i — rx.j\hj) 



(17) 



1 1 <max{/ii , hj } 



Since in this formulation of GSPH there is no need to perform any 
convolution, we implemented Eqs. [T6] and [T7] in the GADGET-3 
code using its original B-spline kernel. Note that the absence of the 
factor y/2 in front of the kernel smoothing length h is due to the 
fact that the above equations have not been derived from the con- 
volution of two Gaussian ker nels, as in the correct GS PH formula- 
tion. FurtheiTnore, following ICha & WhitworthI ( |2003|) . we use for 
this formulation the standard SPH computation of gas density, and 
not the symmetrized one provided by Eq. [15] As we shall demon- 
strate below with hydrodynamical tests, the formulation provided 
by Eqs. l ll6t and Ml\ turns out to be exceedingly diffusive and pro- 
vides an incorrect description of the development of gas-dynamical 
instabilities. This highlights the relevance of properly describing 
the volume convolution in the particle description of the equations 
of fluido-dynamics. 



3.2 The Riemann solver and the slope limiter in GSPH 

We need to know the values of P* and ii* at the interface position, 
for each pair of particles (i, j), in order to use them in the GSPH 
Equations l |13t and il4t . In grid-based hydrodynamical codes, the 
Godunov method is based on solving a Riemann Problem (RP) at 
each cell interface to evaluate numerical fluxes, which are then used 
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to update the value of thermodynamical quantities in the cell. In the 
case of GSPH, we do not need to compute fluxes, since we want to 
preserve the Lagrangian nature of the hydrodynamic description, 
while we only need the values of the (post-shock) velocity v, and 
pressure P, at the interface. 

To solve the RP, we first have to define the right and left (pre- 
shock) states (Pr,l, vr.l and Pr,l)- The simplest choice is to use 
the SPH values of P, v and p, associating the i-th particle to the 
right state and the j-th particle to the left one. This corresponds in 
the original Godunov scheme to a first-order spatial accuracy. Since 
we are not interested in computing fluxes, the exact position of the 
interface in this first-order scheme is not important, given that ther- 
modynamical quantities are assumed to be spatially constant for 
each state. This would make the first-order scheme in principle well 
suited for our purpose. However, this scheme is known to be highly 
dissipative, thus making it not ideal to capture the development of 
fluido-dynamical instabilities. 

A second-order spatial accuracy can be achieved using a 
piecewise linear distribution of thermodynamical quantities. In this 
case, the left and right states are defined by the values of (P, v, 
p) at the interface position s,. The (linear) interpolation of these 
variables is ob tained using their derivatives along the s-axis (see 
llnutsukj|2002i . for details). We describe in the Appendix how the 
position s* of the interface can be computed in the case of linear 
and cubic interpolation of the volume V{s). 

In higher-order (i.e. second-order and above) schemes, ther- 
modynamical quantities must be limited to obtain a stable descrip- 
tion of the discontinuity. This is obtained by implementing a lim- 
iter, which is defined as "a non-linear algorithm that reduces the 
high-derivative cont ent of a subgrid inteipolant i n order to make 
it non-oscillatory" (van Leer 2006). For instance, Inutsuka] 1 I2OO2I) 
imposed that a second-order reconstruction is performed when the 
components of velocity gradients, at the position of the two parti- 
cles, parallel to the s-axis have the same sign, while one resorts to 
a first-order reco nstruction in case of discordant signs. 

We note that Godunov (1959) proved a theorem which states 
that any advection scheme preserving the monotonicity of the solu- 
tion is at most first-order accurate. This holds only if the discretiza- 
tion of the advection equation is linear. Thus, non-linear schemes 
are needed to achieve higher order accuracy. On the other hand, 
for them to be useful, higher-order schemes need to include in the 
interpolator a prescription to limit spurious oscillations. For exam- 
ple, in the context of Eulerian schemes. Van Leer (1979) proposed 
to employ a "harmonic gradient averaging" technique, in which the 
gradient of the thermodynamical quantity Q is the harmonic aver- 
age of the gradients in the k + 1/2 and k + 3/2 cells, namely: 

where A^Q ~ Qk+1/2 ^ Qfc-1/2 and we assumed unity value for 
the cell size. In this notation, Qk+1/2 indicates the mass-averaged 
value of a given thermodynamical variable in the cell k + 1/2, hav- 
ing boundaries at k and + 1, thus with the obvious extension for 
the meaning of Qk-i/2- The mass-averaged gradients of the vari- 
able Q at half-cell positions, A(3fc±i/2, are then used to estimate 
the right-state and left-state values of Q at the interface: 

Qfc±i/2 =F AQfc±i/2/2 if AQfe_i/2AQfe+i/2 > 0, 
Qk±i/2 if AQfc„i/2AQfc+i/2 < 0, 

It can be shown ivan Leej2006l) that such an interpolator cor- 
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Figure 1. Cartoon illustration of how the gradient of a thermodynamical 
quantity Q is computed f or the purpose of implementing the limiter intro- 
duced bv lvan Leei[il97Sj) . The upper figure describes the standard case of 
the computation in case the variable Q is assigned on a regular grid. The 
lower figure is for the case, relevant for our implementation of GSPH, in 
which the vaiiable Q is assigned at the positions of the j-th and j-lh parti- 
cles. 



responds to a standard central difference of the quantities Q, lim- 
ited by a term of order (As)^ which depends on the rate of change 
of the quantity itself through its second-order derivative. Thus, the 
harmonic gradient averaging is a good example of a non linear in- 
terpolator with a built-in limiter. 

In the case of an Eulerian scheme implemented on a regular 
Cartesian grid, Eq. [18] uses three adjacent cells, namely k ~ 1/2, 
fc + 1/2 and fc + 3/2. Clearly, this prescription to compute gradients 
can not be directly generalised to our GSPH scheme, where the 
values of the variable Q are assigned at the positions of the i-th and 
j-th particles, while there is no a third particle to form a regular 
grid along the direction of the separation vector Xi — Xj . 

In the upper part of Figure[T]we schematically show how quan- 
tities are evaluated in the Eulerian scheme, whil e the lower par t 
shows our extension of the implementation of the Ivan Leeild 19791) 
limiter in the case of the unstructured grid, which is relevant for the 
GSPH. In this case, we only have quantities Qi and Qj evaluated 
(19)at the particles positions. We define the equivalent of the Eulerian 
A(5fe_i/2 as: 



AQi = 



As 



(20) 
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where As it the distance between the two particles. To evaluate 
the equivalent of AQj.+i/2, we use a "ghost" particle (j), which 
is located the same distance As from the i-th particle, but in the 
opposite direction along the s axis. Since we have the SPH esti- 
mate of the gradient (VQ)^ at the position of the i-th particle, we 
can compute the expected value Q[j) at the position of the "ghost" 
particle as 



Q{0) ^Q^ + (VQ)z ■ SAS. 



(21) 



Thus, the equivalent of the quantity that we defined in 

the case of the regular grid, is 



AQ2 = 



Qi + (VQ)i ■ s ■ As - 
As 



VQ - s 



(22) 



Therefore, A(32 is simply the projection of the gradient of the 
quantity Q on the s-axis. In general, it is AQi 7^ A(52, since 
Q2 depends on the 3D properties of the field Q through its gradi- 
ent computed at the i-th particle, while the latter only depends on 
the values of the field at the positions of the i-th and j-th particles. 
Of course, the same calculation can be repeated for the j-th parti- 
cle using a "ghost" particle (i), which is defined in the same way. 
Therefore, the extension of Eg. (118) to the case of an unstructured 
grid, defined by particle positions, reads 

Tq11\%\ ifAQiAQ2>0, 



aq 
as ' 



if AQ1AQ2 < 0, 



(23) 



Having computed the values of the gradients ^ at the po- 
sitions of the i-th and j-th particles, we can assign the values of 
thermodynamical quantities on the left and right states at the in- 
terface, following the same procedure as in 102. We will refer to 
this implem entation as "sec ond order reconstruction with van Leer 
limiter" (see lvan Leeil2006L for a detailed description) . 

Once we have thermodynamical variables assigned on the left 
and right states (Pr^l, vr,l, Pr.l), we still have to solv e the 
RP. M any Riemann solvers do exist in the literature (e.g. iTord 
and references t herein). Here, w e use the iterative solver 
originally proposed by ( van Leeill 19791) and described in detail by 



ICha & WhitworfliI tOOJi We refer to this paper for a complete de- 
scription of the implementation of this solver The general idea is 
( i) to define the Lagrangian shock speed W starting from an initial 
guess based on the values of P, and i;,; (ii) calculate the values of 
the tangential slopes Z — ; ( Hi) use these slopes to evaluate the 
new values of P* and v,; (iv) iterate until the variation in P, with 
respect to the previous iteration falls below a given threshold value, 
namely, is less than 1.5 %. We also checked th at using a differ - 
ent solver, namely a Newton iterative solver (e.g. Ivan Leeij|2006h . 
the results of our hydrodynamical tests do not appreciably change. 
Both the Van Leer and the Newton solvers are exact. Usually, they 
provide a converged result in 5-7 iterations. A significant speed-up 
of the code could be obtained by using instead an approximate one- 
iteration solver, such as the Harten-L ax-van Leer-Contact (HLLC) 
solver proposed by dToro et al.lll994 . 

In summaiy, the implementation of the Godunov method in 
our GSPH version of the GADGET-3 code requires the following 
steps. 

1. Estimates of the volume function V{s) and of the position s, 
of the interface between each pair of particles. In our implementa- 
tion, this can be done through either a linear interpolation or a cubic 
spline; 



2. Choice of the spatial order of the reconstruction of the ther- 
modynamical values at the interface position. Currently, we have 
implemented first and second order reconstructions. A third order 
scheme, such as the Piece-wise Parabolic Method (PPM), could in 
principle be implemented. 

3. Choice of the limiter, in the case of second-order reconstruc- 
tio n. We hav e implemented both the "standard" reconstruction, as 
in llnu tsuka (2002), and the reconstruction scheme proposed by 
Ivan Lee r (2006). 

4. Choice of the solver for the RP. We have implemented two 
equivalen t exac t solvers, namely the iterative solver proposed by 
Ivan Leerl d 19791) and the Newton solver. We will present results 
based on the former solver, while we have verified that identical 
results are obtained by using the Newton solver. 



We provide in Table |3.2l a description of the variants of the hy- 
drodynamical schemes that we compare through the tests described 
in the following section. 

The first one is the original GADGET-3 entropy conserving 
scheme, tagged GADGET. We then use a traditional, energy con- 
serving SPH scheme (TradSPH). In this case, however, we use 
a Gaussian kernel, rather than the B-spline one implemented in 
GADGET. We tag as GSPH our new reference implementation, 
which uses cubic spline interpolation for volumes, second-order 
reconstruction of thermodynamical variables at the interface, with 
the limiter a nd the iterative s olver for the Riemann problem, both 
proposed bv Ivan Leerl (Il979l) . Furthermore, we tag as GSPH-I02 
the version based instead on using the limiter adopted in 102, with 
GSPH- 1 ORD the version based on a first-order reconstruction of the 
thermodynamical quantities at the RP interface, and with GSPH- 
VLiN the version based on the linear interpolation for the volume 
function V{s) (see Appendix). Finally, we also implemented th e 
version of the Godunov SPH scheme dCha & WhitworthI |2003|) . 
which is described by Eqs.(|16) and ( |17t , rather than by actual equa- 
tions of GSPH involving the volume integrals, as in Eqs.l[5]l and 
We will refer to this scheme as GSPH-CW. 



4 RESULTS 

We describe in this section the tests of our GSPH implementation 
in the GADGET code. The hydrodynamic tests performed are the 
shock tube, the development of Kelvin-Helmholtz instabilities in a 
shear flow and the disruption of a cold dense blob moving in a hot 
wind. 



4.1 Shock tube 



We first consider the standard Sod shock-tube test ( ISodlll978l) . 
which provides a mean to validate the code capability to describe 
ba sic hydrod)fnaini c features. Initial conditions are the same used 
by ISpringell (l2005l) to test the GADGET SPH scheme. An ideal 
gas with politropic index 7 = 1.4 is considered initially at rest, 
filling half space with gas at unit pressure and density (pi — 1, 
Pi = 1), and the other half space with lower density (p2 ~ 0.25) 
and lower pressure (P2 = 0.1795) gas. Despite the intrinsic one- 
dimensional nature of the test, initial conditions are generated in 
three dimensions with an irregular glass-like distribution of equal- 
mass gas particles. A periodic box was chosen having a longer size 
in the a;-direction, with {Lx, Ly, L^) = (60, 1, 1). A total number 
of 75000 particles have been included in the initial conditions. We 
run this test with four different hydrodynamic schemes (see Table 
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Name 


Hydrodynamic scheme 


Implementation details 


Sod 


KH 


Blob 


GADGET 


SPH 


Entropy conserving with B-spline kernel. 


100 


442 


50 


TradSPH 


SPH 


Energy conserving with Gaussian kernel. 


100 






GSPH 


GSPH 


Based on Eqs. I16t and I17t with cubic spline for volume interpolation, 
second-order reconstruction and limiter bv van Leer (1979). 


100 


100 & 300 


200 


GSPH-I02 


GSPH 


The same as GSPH but based on the limiter bv Inutsuka (2002). 


100 


100 & 300 


200 


GSPH-lORD 


GSPH 


The same as GSPH-I02 but using the first-order reconstruction 
for the RP solution. 


100 


300 


200 


GSPH-VLIN 


GSPH 


The same as GSPH, but using the linear interpolation for the volume 
function V(s), as in Eq. )25t. 




300 


200 


GSPH-CW 


GSPH 


Same as GSPH but based on Eas. fT6l and flTb bv Cha & Whitworth (2003). 


100 




200 



Table 1. Chai'acteristics of the hydrodynamical schemes implemented and of the hydrodynamic tests carried out. Column 1 : name of each scheme; Column 2: 
hydrodynamic method on which each scheme is based; Column 3: basic description of the implementation details of each scheme (see Section|3]for further 
details); Columns 4, 5 and 6: number of neighbour used for each scheme in the Sod shock tube, Kelvin-Helmholtz and blob test, respectively. 



[I2l >. namely GADGET, TradSPH, GSPH, GSPH-CW and GSPH- 
102. In all cases, the test has been mn by using 100 neighbours 
within the kernel. 

We show our results for the SOD test by restricting the range 
of X coordinates to vary between 20 and 40. In fact, since we use 
periodic initial conditions, we have two shocks inside the compu- 
tational domain. We only focus on one of them, while discarding 
the uninteresting unperturbed regions [0:20] and [40:50] along the 
X-axis. Results for this test are shows as scatter plots of density, 
pressure, velocity and entropy, all expressed in internal code units. 
In fact, performing a binning often smooths out details which are 
useful to understand the differences in the behaviours of our the 
different hydrodynamic schemes. 

We show in Figure |2] density, pressure, and velocity at the 
time T = ffl when the shock is well developed, for the GAD- 
GET, TradSPH, GSPH and GSPH-CW schemes (red, green, blue 
and magenta points, respectively). 

A well known problem of SPH codes in solving the Sod shock 
tube is the appearance of a pressure discontinuity ("pressure blip") 
at the position of the contact discontinuity. This is at variance with 
respect to the exact analytic solution, that predict instead pressure 
to be continuous across this discontinuity. This pressure discontinu- 
ity arises as a consequence of the error that the SPH scheme makes 
in correctly estimating the density gradient at the density discon- 
tinuity. As a result, a sort of spurious surface tension force arises, 
due to the opposite signs that the pressure gradients have on the two 
sides of the discontinuity. In fact, GADGET and TradSPH formu- 
lations clearly show a pressure blip at x ~ 33.5, where the discon- 
tinuity is located. This is emphasised in the bottom left panel of Fig. 
|2l which provides a zoom of the pressure in this region. Since this 
spurious pressure feature arises from errors in the computation of 
density at the discontinuity, its origin can be traced back to the in- 
correct estimate of the 1/p SPH volume associate to each particle. 
Two conceptually different solutions to this problem can be then 
devised. A first one is based on improving the density estimate in 
the computation of density gradients across the discontinuity (e.g. 
iRead et alj|2olol) . A second one relies instead on the introduction 
of thermal diffusion across the discontinuity, which masquerade the 
effect of the error in the volu me estimate and prevents the onset of 
the surface tension there (e.g. |Pricel2008l) . 

^ Note that time is a dimensionless quantity in this test, thus it does not 
depend on our chosen systems of units 



ICha & WhitworthI jlOOih noticed that the GSPH-CW scheme 
was in fact able to prevent the development of the density blip. 
Since this scheme does not pay attention to correctly estimate vol- 
umes, its behaviour should be ascribed to the inclusion of a diffu- 
sion term. This term is indeed provided by the [v* — x*] appear- 
ing on the r.h.s. of Ea.J16t. which in fact described a net exchange 
of thermal energy, with a zero total mass exchange, between each 
particle pair. Its effect is similar to that of adding an artificial ther- 
mal diffusion term in the energy equation. The main difference, 
however, is that here such a mixing term naturally arises from the 
Godunov scheme. 

As shown in the bottom-right panel of Fig.|2l we do confirm 
that this scheme does not produce a pressure discontinuity, which 
is replaced by an oscillation around the exact solution. Therefore, 
while diffusion prevents the "pressure blip", inaccuracy in the den- 
sity estimate is still present and induces the appearance of a spuri- 
ous, although much reduced, pressure force at the contact disconti- 
nuity. Quite interestingly, this spurious force is further, and greatly, 
reduced in our reference GSPH scheme. In this case, second-order 
accuracy in density estimate is enforced through the computation 
of the convolution integrals in Eqs. (|5j and ([8} (see also 102). The 
net effect is a significant improvement in the behaviour of pressure 
across the discontinuity. 

To further emphasise the different behaviour of the SPH-based 
(GADGET and TradSPH) and GSPH-based (GSPH and GSPH- 
CW) schemes, we plot in Figure |3] the fractional variation of the 
entropy of particles as a function of their final position, across the 
discontinuity. Quite clearly, entropy variation is only determined 
by the presence of weak shocks for the GADGET and TradSPH 
scheme and, therefore, can only have positive sign. On the other 
hand, the presence of diffusion in the GSPH and GSPH-CW allows 
an exchange of internal energy across the discontinuity. This dif- 
fusion manifests itself with both positive and negative variations of 
entropy, with lower-density particles located on the right side of the 
discontinuity loosing thermal energy in favour of particles located 
on the other side of the discontinuity. It is such a mixing that is 
responsible of the damping of the pressure blip. 

Both GADGET and TradSPH schemes show pressure wiggles 
immediately before the shock wave, which is located at 35 < a; < 
37 (see the inset in the pressure panel in Fig.|2). Such wiggles in 
pressure correspond to wiggles in the velocity, as shown in the top- 
right panel of Fig|2] at the same position. Note that velocity wig- 
gles are largest for GADGET and smallest for GSPH. Wiggles are 
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20 25 30 35 40 25 30 35 40 




20 25 30 35 40 32.5 33 33.5 34 34.5 35 



Figure 2. Comparison between the results of different hydrodynamic schemes for the Sod shock tube. Scatter plots of the particle density, velocity and 
pressure are shown in the upper left, upper right and bottom left panels, respectively. The bottom right panel shows a zoom of the pressure around the contact 
discontinuity. Red points are for the entropy-conserving GADGET fomiulation, green points for the TradSPH formulation, blue points for the GSPH-CW 
formulation and magenta points for our reference GSPH formulation (see Table |3.2V In each panel, the black continuous line indicates the exact solution. The 
insets emphasise the different behaviour of the four schemes, as also discussed in the text. 



also present in the GADGET run at the same position for the den- 
sity variable (upper-left panel). The insets in the two panels show a 
blow-up of the region around the shock wave positions, and empha- 
sise the presence of such wiggles. The prominence of the pressure 
wiggles in the GADGET simulation is due to the lack of thermal 
energy diffusion that characterises this hydrodynamic scheme. In a 
sense, these wiggles have the same origin as the wiggles in veloc- 
ity appearing i n the shock-tube test when artificial viscosity is not 
included (e.g.. lRosswog..2009.) . In fact, artificial viscosity removes 



velocity wiggles since it effectively provides a diffusion of momen- 
tum at the shocks. In a similar way, diffusion of thermal energy as- 
sociated to the GSPH scheme is effective in removing wiggles in 
pressure at such discontinuities. 

In the upper-right panel of Figure |2] we also notice that the 
scatter in the velocities is minimum for GSPH and maximum for 
GADGET. The shock wave is better captured by GADGET and 
GSPH, with GSPH-CW performing worse. It is quite remarkable 
that GSPH is able to correctly capture the shock, while preventing 
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Figure 3. Entropy variation in GADGET, TradSPH, GSPH and GSPH- 
CW hydrodynamic scheme for the Sod shock test around the contact dis- 
continuity. Colour coding is the same as in Figure |2] Here we show the 
difference between the final and initial entropy of particles, normalised to 
their initial entropy, as a function of final particle positions. 



at the same time the appearance of noise in the velocity field, with- 
out the introduction of an artificial viscosity term in the momentum 
equation. Furthermore, the better accuracy of GSPH with respect 
to GSPH-CW is the consequence of the improved accuracy of the 
former scheme. 

The shape of the rarefaction fan in density and pressure (up- 
per and lower left panels in Fig. 1, respectively) are better captured 
by GADGET and TradSPH formulations, while GSPH-CW shows 
a much smoother behaviour. GSPH also shows a slight smoothing 
of density and pressure at the beginning and at the end of the rar- 
efaction fan. On the other hand, GADGET produces a piling-up of 
particles at the onset of the fan, corresponding to an excess of par- 
ticles having low velocities, at the same position. This feature can 
be better appreciated in the lower inset of upper-left panel of Fig. 
1. 

Figure |4] shows results for the Sod test, when we use different 
implementation for the GSPH formulation. Besides our standard 
implementation of GSPH, we also show results based on the lim- 
iter used by 102 in its implementation of GSPH (GSPH-I02) and on 
the first-order reconstruction at the interface s, for the solution of 
the RP (GSPH-Iord; see Tabl d3.2b . From the behaviour of density 
and pressure, it is clear that a first-order reconstruction results in a 
more diffusive behaviour: the shape of the rarefaction ramp and of 
the velocity profile at the position of the shock wave are smoother 
than for the other two implementations. For the same reason, no 
velocity wiggles appears when this formulation is used. The dif- 
ference between the two limiters is instead clear when we anal- 
yse the density right before the rarefaction fan. There, the standard 
limiter produces an accumulation of particles, similar to what we 
found when using the SPH GADGET formulation. Such an accu- 
mulati on is instead eliminated by adopting the limiter bv Ivan Leer! 
( Il979h . Furthermore, the use of this limiter avoids the accumula- 
tion of low-velocity particles after the position of the shock wave. 



as can be seen in the upper-right panel. The pressure blip at the 
density discontinuity is erased in all of these three schemes, thus 
we do not show a zoom-in for the pressure. 



4.2 Kelvin-Helmholtz instability 

The Kelvin-Helmholtz instability occurs across a contact disconti- 
nuity in the presence of a tangential shear flow. It belongs to a class 
of tests that have been used over the last few years to assess the 
capability of SPH and grid-based methods to capture fluidodyneim- 
ical instabilities (e.g. Agertz et al. 2007; Price 2008; wadslev et aU 
[2008 ; Read et al...201Q ; .Valcke et al.,„201Q) . It describes the evolu- 
tion of two fluids having different densities and in pressure equi- 
librium, moving with opposing velocities. The interface between 
the fluids is perturbed leading to a phase in which fluid layers de- 
velop vortex instabiliti es, with subsequent mixing of the two phases 
jChandrasekhad [T96II) . We remind here that a two-dimensional 
version of this t est has been also recently discussed in detail by 
IChaetalJjioiol) to assess the capability of GSPH to describe insta- 
bilities. We extend here this test in three-dimensions for our GAD- 
GET implementation of GSPH. 

The KH test that we present here belongs to the "Wen- 
ge n" suite of hydrod ynamical testfl and is described in detail 
by iRead et al.1 ilOUt) . An ideal idea gas with politropic index 
7 = 1.4 and mean molecular weight /i = 1 is assumed. Ini- 
tial conditions are generated in the periodic simulation domain 
{Lx, Ly, Lz) — (256,256, 16) kpc centered on the origin. Two 
domains with \y\ < 64 kpc and \y\ > 64 kpc corresponds to the 
two fluids having pi = 2p2, Ti = O.5T2 and opposing velocities 
with the same modulus v. In this test, p2 = 3.13- 10~* M0 kpc""^, 
T2 = 3 ■ lO^K, and « = 40 km s~^. Equal-mass gas particles are 
initially located on a grid, whose spacing in the two domains is set 
in such a way to reproduce the difference in density. Instabilities are 
triggered by imposing a velocity perturbation along the j/-direction, 
having a characteristic wavelength A = 128 kpc. The characteristic 
KH time-scale for the development of instabilities from this pertur- 
bation is 



2v{pip2Y'^ 



(24) 



The units used in the Wengen tests are: kpc for length, km s^^ 
for velocities, and 10^° for masses. In this system, the unit of 
time is t, = 0.977 Gyr. Using the perturbation described above, 
tkh — 3.32 Gyr. Initial conditions have been generated using 
774. 144 gas particles. We carried out this test using different im- 
plementations of the GSPH scheme: our reference scheme (GSPH), 
the scheme based on the limiter adopted by 102 (GSPH-I02), the 
scheme based on a first-order reconstruction for the thermodynam- 
ical quantities at the interface (GSPH-Iord), and the scheme based 
on the linear interpolation of the volume function V{s) (GSPH- 
VLin). In order to reproduce the results reported for the GADGET 
code on the web site of the Wengen tests, we carried out the KH 
test in this case using 442 neighbours. As for the different imple- 
mentations of the GSPH scheme, we always used 300 neighbours 
within the Gaussian kernel, while we also checked the effect of us- 
ing instead 100 neighbours for the reference GSPH scheme and for 
the GSPH-I02 scheme. 

We show in Figure |5] the development of the KH instabil- 
ity for GADGET, GSPH and GSPH-I02at three different times. 



|http : / /www ■ astrosim. net/ code| 
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Figure 4. Compaiison between the results of different implementations of the GSPH Equations (5j and (Sj for the Sod shock tube. Scatter plots of particle 
density, velocity and pressure are shown in the upper left, upper right and bottom panels, respectively. T he bottom - right panel shows a zoom-in of pressure 
behaviour around the onset of the rarefaction fan. Red points are for the GSPH with the standard limiter of llnutsukal ( |2Q02|) (GSPH-I02), green points for the 
same scheme when using a first-order reconstraction of density and pressure at the interface in t he solution of the Riemann problem (GSPH-lORD), while 
the blue points correspond to our standard implementation (GSPH) based on using the limiter by Ivan Lee3 j 19971) and the third-order spline interpolation of 
volumes. 



namely 0.5tkh,tkh and 2tkh- It is clear that, while the stan- 
dard entropy-conserving SPH scheme dumps the instability, both 
our GSPH and GSPH-I02 schemes successfully capture its devel- 
opment. Alt = tkh, vortexes begin to show up, and alt — 2tkH 
they are fully developed and display the typical "cat-eye" structures 
in gas density. 



This figure confirms the results reported by ICha et all (1201 Ol) 
and demonstrates the capability of the GSPH scheme to develop 
gas-dynamical instabilities. Note that, while the damping of the 



instability in the GADGET entropy conserving scheme is due to 
the use of an artificial viscosity and to development of a "artifi- 
cial surface tension" force at the interface, originat ed by the re - 
pulsion of SPH particle at the discontinuity (see e.g lPricd ( l2008h : 
ICha et all (jlOlCf)), the former is absent and the latter strongly re- 
duced in GSPH schemes (see Inutsuka (2002); Chaetal. (20i3) 
for a discussion of the different estimate of density on GSPH and 
on its effect on the artificial surface tension). The much improved 
description that the Godunov SPH scheme provides in describing 
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Figure 7. Pressure of gas particles as a function of tlieir y coordinates for 
GADGET (red dots) and GSPH (green dots), at time t = 2tkh, in cor- 
respondence of tlie same vortex sliown in Fig. |6] 78 < x < 210 kpc, 
161< y < 237 kpc, -3.5 < 2 < 5.5 kpc. 



the development of KH instabilities has to be ascribed to the fact 
that this method is based on explicitly computing the convolution 
integrals in Eqs. ([5} and ([8}, to a 0{h?) accuracy, through the in- 
terpolation of the volume function V{s). 

Figure |6] shows a blow-up of one of the curl structures, which 
forms after the KH instability is fully developed, for both the GSPH 
and the GSPH-I02 schemes. The difference is small, but the Van 
Leer limiter is able to more neatly capture the vortex structure. We 
hypotize that the reason for this behaviour is that the Van Leer lim- 
iter is capable to slightly reduce the intrinsic residual numerical dif- 
fusion associated to the Riemann solver, with respect to the limiter 
implemented by 102. 

To further highlight the different ways in which entropy- 
conserving SPH and GSPH respond to the velocity perturbation 
across the contact discontinuity in the KH test, we show in Fig- 
ure [7] a scatter plot of the pressure of gas particles as a function 
of the y coordinate (i.e. in the direction parallel to the direction to 
the velocity perturbation triggering the KH instability). This scatter 
plot includes only the particles contained within the same vortex 
region shown in Fig. |6] It is clear that GADGET develops a pres- 
sure "blip" at the discontinuity, similarly to what happens in the 
shock tube test. On the contrary, the blip basically disappeared in 
the GSPH scheme, and a complex pressure structure, associated to 
the vortex, develops. 

In order to further highlight the role that different details in 
the implementation of the GSPH scheme have in the development 
of the KH instability, we also studied the effect of varying the 
number of neighbours, the reconstruction prescription and the vol- 
ume estimates. The results of these tests are shown in Figure [8] 
The most striking effect is given by the reconstruction order for 
the thermodynamical quantities at the interface where we solve the 
Riemann problem. A first-order reconstruction ( GSPH- Iord; upper 
right panel) is so diffusive that the KH instability does not develop 
at all. This demonstrates how crucial it is for the development of the 



KH instability to choose a scheme that gives the lowest possible de- 
gree of diffusion where it is not required. In fact, the GSPH- Iord 
has been shown to be effective in preventing the formation of the 
pressure blip in the shock tube test (see Fig|2). However, the ex- 
cess of diffusion, which manifests itself in the shock tube test as a 
smooth transition at the rarefaction fan, is such to prevent the devel- 
opment of the KH instability. We also note that it is quite important 
to use a rather large number of neighbours: decreasing its number 
from 300 to 100 (upper central panel) produces a smoother struc- 
ture of the vortexes, and fails to capture their "cat-eye" shape. This 
is due to the increase of noise in the density estimate at the discon- 
tinuity, that arises when a smaller number of neighbours is used. 
On the other hand, the KH test is also sensitive to the precision of 
interpolation of the volume function: using a linear interpolation of 
this function (GSPH-Iord; lower central panel), develops instabili- 
ties whose "curl" structure is however less resolved than for a cubic 
interpolation. As for the GSPH-CWscheme, based on Eqs. ( I16t and 

the large amount of diffusivity is such to prevent the devel- 
opment of KH instabilities. The lack of accuracy of this scheme is 
due to two main differences with respect to our reference GSPH 
scheme: firstly, Eqs. ( 116) and ( 117) are not derived from the convo- 
lution of the equation of motion and energy conservation; secondly, 
in this scheme we have no means to locate the position of the inter- 
face for the solution of the RP, thus implying that we are effectively 
resorting to a first-order reconstruction. 

As a further check, we run another KH test, in which we mul- 
tiplied the y-axis velocity perturbation by a factor of five. The aim 
of this test is to verify whether lower-order schemes succeed in 
following the development of the KH instability when the pertur- 
bation is strong enough. We show the results of this test in Figure|9] 
at the time t = tkh of our "standard" test, so that also the ampli- 
fication of the perturbation is appreciable. In this case, the entropy 
conserving GADGET scheme develops arms, but fails to capture 
the development of the vortexes. On the other hand, GSPH is con- 
firmed to successfully capture the instability, while the GSPH-Iord 
and GSPH-CW schemes confirm to be very diffusive and to smooth 
out the instability. The results of this test demonstrate that a higher- 
order GSPH scheme is necessary to correctly treat the development 
of KH instabilities. 

In summary, the results shown i n this Section co nfirm and ex- 
tend the 2D test results presented bv lCha et al.l ( |20IO|) on the capa- 
bility of GSPH to follow the development of KH instabilities. We 
demonstrated that the per formance of GS PH is further improved 
by adopting the limiter bv Ivan Leer! ( Il979h . instead of that of 102, 
used by Cha et al. ( 2(513). Furthermore, our results also highlight 
that the development of the instability is inhibited in different ways 
hy (a) errors in density estimate at the discontinuity, as in standard 
SPH or when using a small number of neighbours in GSPH, and 
(b) numerical diffusion, which increases when using a less accu- 
rate first-order reconstruction of thermodynamical variables at the 
interface. 



4.3 The blob test 

This test describes the disruption of a cold gas cloud having uni- 
form density, moving in pressure equili brium against ajiot lower- 
density wind. This test has been used by A gertz et al. I( l2007l) to as- 
sess the capability of different hydrodynamic schemes to describe 
the blob disruption d ue to the onset o f KH and Rayl eigh-Taylor 
instabilities (see also iRead et al.ll201oh . ICha et al.l mm) recently 
used a two-dimensional version of this test to assess the perfor- 
mance of their GSPH implementation. The version of the blob test 
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Figure 5. Development of the Kelvin-Helmholtz instabilities at different times, as described by the GADGET (upper panels), GSPH-I02 (central panels) and 
GSPH (bottom panels) scheme. Results are shown at three different times, G.bTKH,TKH and 2tkh in the left, middle and right panels, respectively. Each 
panel report the the projected density for a slice of coordinates —3.5 < z < 5.5 kpc. Colour scale is linear and ranges from 2 x 10~^MQkpc~^ (black) to 
5 X lO-''M0kpc-2 (white). 




Figure 6. Blow-up of the curl structure developed by the KH instability at t = 6.75, for GSPH (left panel) and by GSPH-I02 (right panel). The structure 
shown here developed in the ranges of coordinates 78 < a; < 210 kpc and 161 < y < 237 kpc. We show the projected density for a sUce of coordinates 
—3.5 < z < 5.5 kpc. Colour scale is linear and ranges from 2 x lO~''^M0kpc~^ (black) to 5 x lO~'''M0kpc~^ (white). 
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Figure 8. Effect of different GSPH implementations on the KH instability at t = 2tkh- Results are shown for our reference GSPH scheme using 300 
neighbours (upper left panel) and 100 neighbours (upper central panel), using GSPH- lORD with first-order reconstruction at the interface (upper right panel), 
and using GSPH-I02 with the limited by 102 (bottom left panel), using GSPH-VLIN with linear interpolation for the volume function (bottom central panel), 
and for the GSPH-CW scheme based on Eqs. (TTJ and )16t (bottom right panel). Ranges of coordinates and gray-scale density coding are the same as in Fig. 

13 



that we use here also belongs to the "Wenge n" suite. A full de - 
sciiption of the initial conditions is provided bv lReadetallbOld) . 
The simulation domain is given by a periodic rectangular box with 
{Lx, Ly, Lz) = (2, 2, 8) Mpc, with the origin of the coordinates 
located at the centre of this domain and the blob initially located 
at (0, 0, —3) Mpc. The radius of the cloud is set to Vc — 197 kpc. 
Internal density of the cloud is a factor x = 10 higher that in the ex- 
ternal medium, with the temperature being correspondingly a factor 
10 lower so as to fulfil the condition of pressure equilibrium. Initial 
conditions are generated by placing equal-mass particles in a lattice 
configuration, so as to satisfy the above density requirements. The 
velocity of the wind is = 1000 kms^^. An initial instability is 
also used to the surface layer of the cloud to trigger a large-scale 
instability (see Read et al. 2010, for a full description of the initial 
conditions). Units of measure for length, velocity, mass and time 
are the same as in the KH test. Initial conditions are generated using 
10^ gas particles. We use 200 neighbours for the test based on the 
different implementations of GSPH, and 50 for GA DGET, which 
uses the B-spline kernel (see Table 13. 2t . Following lAgertz et al.l 
( l2007h . we define the cloud crushing time, Tcc = = 0.61 

Gyr, which gives the typical time-scale for the evolution of the 
cloud moving at supersonic velocity. 

The results of this test have implications for a number of rele- 
vant astrophysical and cosmological applications, in which a dense 
gas cloud interact with a lower density medium. For instance, this 
is the case of a cold molecular cloud in the inter-stellar medium, 
which interacts with ejecta from a nearby exploding supernova. 
Another example is provided in cosmological simulations by sub- 
structures bringing relatively cold gas which merge into larger ha- 



los permeated by hotter gas during the hierarchical assembly of 
galaxy groups and clusters. 

The blob is expected to be initially destabilised by Richtmyer- 
Meshkov instability (RM) and by Rayleigh-Taylor (RT) instabil- 
ity and subsequently further dissolved by KH instability. We show 
in Figure [TT]the evolution of the projected gas density in this test 
at three different times, for the GADGET (upper panels), GSPH 
(bottom panels) and GSPH-I02 (middle panels) hydrodynamical 
schemes. As expected, the limitations of the standard entropy con- 
serving formulation of SPH in following hydrodynamical instabili- 
ties make this scheme unable to describe the disruption of the cold 
blob (see lAgertz et al . 2007, for a detailed discussion). Already at 
early times, t — 4 (left panels), the GADGET simulation develops 
less hydrodynamical instabilities in the up-wind part of the blob. 

Even if the total number of gas particles in our blob simula- 
tion is fairly high, force resolution, which is directly related to the 
mean interparticle separation, is lower than in the KH simulations. 
In fact, the mean interparticle separation is 14.7 kpc for the blob 
test, and 0.89 kpc for the KH test. For this reason, the development 
of hydrodynamical instabilities is not clearly visible in Figure [TT] 
In order to show the different behaviour of the different numerical 
schemes, we show in Figure [Tol the velocity field in a thin slice, 
centred on the blob, at the time t — 4 when such instabilities begin 
to develop. The upper-left panel shows the field for GADGET, the 
upper-right one for GSPH , the lower-left panel for GSPH-CW and 
the lower-right one for GSPH-Iord. Velocities are computed with 
respect to the rest-frame of the blob. 

The difference between GADGET and GSPH is striking. Hy- 
drodynamical instabilities create vortexes in the back of the blob 
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Figure 9. KH test with a stronger velocity perturbation on the y axis. We show the result at the time tkh = 1 of our standard KH test. Results are shown 
for the entropy conserving GADGET scheme (upper right panel), for our reference one GSPH (upper left panel), for GSPH-CW (lower left panel) and for 
GSPH- 1 ORD (lower right panel). Ranges of coordinates and gray-scale density coding are the same as in Fig.|5] 



in both cases. Quite clearly, the GSPH scheme is by far most ef- 
fective in resolving these vortexes and creating a downstream ve- 
locity flux. Also, the disruption of the blob front due to RT insta- 
bilities is marginally apparent in this figure. Note that the GSPH- 
CW scheme is also more effective than GADGET in capturing the 
vortex structures. This is due to the absence of an artificial viscos- 
ity in such a scheme, which prevents the development of of vor- 
tex structure in the velocity field. However, the large diffusivity of 
GSPH-CW causes the suppression of the downstream flux. Finally, 
the behaviour of the GSPH-Iord scheme is intermediate between 
the other two GSPH schemes, since it is more diffusive than stan- 
dard GSPH but less so than GSPH-CW . 

The difference of the GADGET evolution with respect to the 
two GSPH implementations shown in Fig.[TT]is even more apparent 
at t = 6. At this time the RT and RM instabilities appearing in the 
GPSH simulations are further dissolved into filamentary and curl- 
like structures, which are produced by the onset of KH instabili- 
ties. At t = 8 (right panels) the bulk of the gas initially contained 
in the cold blob still forms a compact structure in the GADGET 
simulation. On the contrary, in both GSPH simulations the blob is 



basically dissolved at this time. The result for the GSPH-I02 case 
ar e fully consistent with the two-dimensional blob test presented 
bv lChaetalJ ( l2010l) . who also used the 102 limiter. 

A compariso n between the results obtained from the 102 and 
the lvan LeeJ ( 1 1979;) limiters show that they perform quite similarly 
in this test. This suggests that differences between these two lim- 
iters only becomes evident in higher resolution tests, such as the 
KH test shown above. 

In order to further verify the performances on the blob test 
of different implementations of the GSPH scheme, we compare in 
Figure [T3l the results obtained at f = 8 for the reference GSPH 
scheme (upper left panel), for the GSPH-VLin version based on the 
linear interpolation of the volume function (upper right panel), for 
the GSPH- 1 ORD version based on the first-order reconstruction for 
the solution of the RP (bottom left panel), and on the GSPH-CW 
scheme based on Eqs. ( 1161 ) and The result for GSPH-VLin is 
rather similar to that of GSPH, although the latter develops a lower 
degree of filamentary structures and instabilities. This result agrees 
with what shown in Fig. [8] for the KH test. On the other hand, the 
results dramatically change if we use instead the first-order recon- 
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Figure 10. Velocity field, at the time t = 4, for the "blob" test, in the region surrounding the blob itself. We plot a slice near the center of the blob, with 
900 < x < 100 kpc, 400 < J/ < 1600 kpc and 1500 < z < 3500 kpc. We show the velocity field for GADGET(upper-left panel), GSPH(upper- 
right panel), GSPH-CW (lower-left panel) and GSPH-lORD (lower-right panel). Each arrow shows the velocity of one gas particle; we undersampled the 
simulation by a factor 0.0025, for clarity. Velocities are computed in rest frame of the blob. 



struction scheme of GSPH-Iord for the assignment of the thermo- 
dynamical variables at the interface. In line with the KH test result, 
the higher degree of diffusivity of this scheme dumps the develop- 
ment of instabilities, thus preserving the structure of the blob. This 
result highlights that, while using an accurate scheme of interpola- 
tion for the volume function has a sizable effect, a much more dra- 
matic change is provided by properly choosing the reconstruction 
scheme for the RP solution. Finally, the blob test confirms the result 
based on the KH test on the incorrect description of the GSPH-CW 
scheme in describing the development of instabilities. 

In order to better quantify the different efficiency that differ- 
ent schemes have in describing the disruption of the blob, we show 
in Figure [T2l the ev olution of the blob mass loss for our various 
schemes. Following lAgertz et al.l ( l2007h . we define a gas particle 
to belong to the blob whenever its density is p > 0.64pci, being 
pel the blob density as set in the initial conditions. As expected. 



GADGET show the lowe st degree of mass los s, a r esult which is in 
line w ith what shown by lAgertz et al.l ( l2007h and iHeB & Springell 
( I2OIOI) . GSPH , GSPH-VLiN and GSPH-I02 have very similar mass 
loss rates, while GSPH-CW retains more mass at the end of our 
simulation. This can be understood in terms of a lesser ability of 
the latter scheme to capture the development of hydrodynamical 
instabilities, as shown in Figure[TO] We also note that the strongest 
mass loss rate takes place for GSPH- 1 ord. Owing to the quite poor 
peiformance of this scheme in describing the development of the 
KH instability, we argue that numerical diffusion is in this case the 
main driver of the blob mass loss. This is also indicated by the dif- 
ferent final shape of the blob (Figure \13) . The cloud is disrupted 
into several pieces in the GSPH and GSPH-VLin schemes, while it 
retains its integrity in the GSPH-Iord and GSPH-CW ones. There- 
fore, significant mass loss of the blob can be obtained both through 
spurious numerical diffusion and through the development of gen- 
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Figure 12. Mass loss of the blob as a function of time for the different 
hydrodynamical schemes. We define a gas particle to be part of the blob 
if its density is p > 0.64/9(.;. with p^i the initial density of the blob. Red 
line with crosses refer to GSPH, green ones to GSPH-102, blue ones to 
GADGET, magenta ones to GSPH-CW, cyan ones to GSPH-VLIN and 
yellow ones to GSPH-lORD. 

uine instabilities. However, only the latter are capable also to pro- 
duce the disruption of the blob into several pieces. 



5 CONCLUSIONS 

In this paper we presented results of 3D standard hydrodynamical 
tests for different implementations of the Godunov Smoot hed Parti- 
cle H ydrodynamics (GSPH) within the GADGET-3 code dSpringell 
1200^ . The conceptual difference between GSPH scheme and stan- 
dard SPH scheme lies in the fact that the former is based on explic- 
itly convolving momentum and energy equations with the interpo- 



lation kernel. The resulting equations implemented in the simula- 
tion code are exact to 0{h^). Suitable expression for momentum 
and energy equations to be implemented numerically are obtained 
by assuming a Gaussian shape for the interpolating kernel. A nat- 
ural way of implementing the equations of GSPH is b y solving 
the R iemann Problem between each pair of particles (see Inutsuka] 
12001 for a detailed discussion of GSPH). Quite remarkably, solv- 
ing the RP between each particle pair brings the extra benefit that 
no artificial viscosity is required to capture shocks, unlike in the 
standard SPH. 

The different implementations of the GSPH scheme, that we 
presented in this paper (see Table 1), correspond to (a) using either 
linear or cubic-spline interpolation for the volume function, which 
provides the position of the interface where to solve the RP; (b) 
using either a first-order or a second-order reconstruction scheme 
to assign thermodynamical variables at the interface in the solution 
of the RP; (c) using different limiters to prevent oscillations of in- 
terpolated quantities in the RP. Furthermore, we also considered a 
variant of the GSPH scheme, which is not based on the convolu- 
tion of the equations of fluido-dynamics, and is instead essentially 
based on replacing pressure and velocity in the SPH equations with 
the values obtained from the solution of the RP. The performances 
of these different implementations to describe discontinuities and 
development of gas-d ynamical instabilities have been assessed us- 
ing a shock tube test (lSodll978t) , a shear-flow test to follow Kelvin- 
Helmholtz (KH) instabilities and the disruption o f a cold blob mov- 
ing in a hot atmosphere (e.g..lA gertz et al.l2007h . 

The results of our simulation tests can be sutimiarised as fol- 
lows. 

(1) As for the shock tube (see Fig.|2j, we verified that GSPH is 
able to correctly follow the development of the shock, despite the 
fact that it does not include artificial viscosity. Furthermore, GSPH 
is also effective in removing the spurious "pressure blip" generated 
by standard SPH at the contact discontinuity, thanks to its capa- 
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bility to describe diffusion of entropy across the discontinuity (see 
Fig.O. Quite interestingly, tfie best description of the discontinu- 
ity is provided by our reference GSPH scheme, rather than by the 
more diffusive GSPH-CW scheme. This result highlights that in- 
cluding thermal diffusion is not enough in itself to provide a com- 
pletely correct description of discontinuities. Indeed, the more ac- 
curate description of density gradients offered by the GSPH, also 
significantly contribute to suppress spurious pressure forces at such 
discontinuities. 

(2) Unlike standard SPH, our reference GSPH scheme is quite 
effective in following the development of KH instabilities in the 
shear-flow test (see Fig. |5}. We have verified that the accuracy in 
developing the curl structure of the instability is quite sensitive to 
the details of the implementation (see Fig. [8). For instance, using 
a first-order reconstruction dramatically degrades the GSPH per- 
formance for this test. Also using a too small number of neigh- 
bo urs, a linear inte rpolation of the vol ume func t ion an d the limiter 
bv llnutsukal ( l2002h . instead of that bv lvan Lee3 ( Il979l) . also some- 
what worsen, although to different degrees, the description of the 
KH instabilities. 

(3) Similar results also hold for the "blob test". Also in this case, 
our standard GSPH implementation follows the onset of Rayleigh- 
Taylor and Richtmyer-Meshkov (RT, RM) and KH instabilities. 
As a result, the blob is dissolved much more efficiently than in 
SPH simulations. However, the performance of GSPH significantly 
worsen in case a first-order reconstruction scheme is used to assign 
variables at the interface. This highlights once again that diffusivity 
of the solution of the RP needs to be minimised to reliably follow 
the development of instabilities in GSPH. 

In this paper we focused on the comparison between different 
hydrodynamical schemes, when applied to control test cases, with- 
out analysing the behaviour of each of such schemes when reso- 
lution is progressively increased (or degraded). On the other hand, 
it is worth reminding that a numerically diffusive scheme could 
converge to the correct solution when ap plied to test cases with 
sufficiently high resolution. For instance, iRobertson et al.l ( 1201 Ol) 
showed this to be case for the Galilean invariance in Eulerian codes, 
in the case of KH instabilities. 

In general our results agree with and extend those presented 



bv lChaetalJ ( l201(]() on two-dimensional KH and "blob" test of the 
GSPH. Our analysis highlights the important role played by recon- 
struction at the interface, by the choice of the limiter and by the 
interpolation order for the volume function (see Appendix). The 
remarkable improvement shown by GSPH with respect to the stan- 
dard SPH to describe contact discontinuities and development of 
gas-dynamical instabilities makes in principle this hydrodynamic 
scheme highly promising for applications in computational astro- 
physics and cosmology. 

As for its computational cost, it is worth pointing out that the 
need of solving the RP between each pair of particles does not rep- 
resent a limiting factor. Clearly, the solution of the RP requires an 
iterative procedure, which in principle could increase the computa- 
tional cost. This could be avoided using an approximate Riemann 
solver, which does not require an iterative procedure, as e.g. an 
Harten-Lax-van Leer-Contact (HLLC) solver. 

However, the lack of artificial viscosity in the GSPH makes the 
Courant condition much less stringent in the shock regions than for 
SPH, thereby leading to a more relaxed time-stepping. We verified 
in our test that this compensates the overhead associated to the RP 
solution. 

Ho wever, it i s wort h remind ing that the GSP H equations de- 
rived bv llnutsukal ( l2002h (see also lCha et ai]|2010 l) and used in our 
implementation (Eas.ll3landll4t hold only for a Gaussian kernel. 
The subsequent request for a fairly large number of neighbours and 
the \/2 multiplicative factor in front of the kernel smoothing length 
in the above equations make the neighbour search quite expensive. 
An implementation of GSPH based instead on a kernel with com- 
pact support, like the B-spline kernel, would clearly be highly de- 
sirable to make the code more efficient for applications involving 
large dynamic and temporal ranges. 
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APPENDIX. VOLUME INTERPOLATION 

In this Appendix we summarise for completeness the expressions 
for the inteipolating volume Vi^j(h), which appears in the GSPH 
equations of evolution JlOt and Jilt , in the case of linear and of 
cubic spline interpolation. We will also provide the corresponding 
expressions for the position of the interface at which the Riemann 
problem between the i-th and the j-th particle is solved. A full 
derivation of all such expressions is provided by 102. 

The linearly-interpolated expression of the specific volume at 
the coordinate s, along the axis joining the i-th and the j-th particle, 
is 



V{a) ^ p{s)-^ =C,,ja + Di,, 
where 



(25) 



A, 



\/(xo + nxi) 



(26) 



We remind that we denote with Si and Sj the components of the 
and Xj vectors along the s-axis, so that Asij = Si — Sj — jx^ — Xj|. 

Including the above expression for p^^ into the integral ap- 
pearmg on the l.h.s. of Eo.(fT2t. one obtains 



V^^'j ^ ^h^Cl, + Dl, . ill) 

In order to compute the position of the interface, let us define the 
weighted-average f*j of a generic function / (x) through the rela- 
tion 



/(x) 

p2(x) 



/ i,3 



W(x - X,; h)W{yi - Xj ; /i)dx = 



■W{yL - X,; h)W{yi - Xj ; /i)dx . 



p2(x. 



(28) 



Using then the linear approximation for /(s) = a{fi — fj)/Asij 
and the above linear interpolation for p(s)~^, one obtains 



r* _ fi fj * I /' + /j 



(29) 



(30) 



In the above equation the position of the interface 

is defined as the position on the s-axis at which the linearly- 
interpolated function / takes the value f*j . 

The computation for a more accurate cubic-spline interpola- 
tion of the volumes proceeds in a similar way. In this case, it is 



V{s) ^ p ^{s) = Aijs +BijS +Cijs + Di,j 



(31) 



Eqs.(61) of 102. Using again Eq.lll2ll for the definition of Vij, one 
obtains in this case 



+ ^h\2B,,,D,.,,+Cl) + Dl 
so that the expression for the position of the interface becomes 
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(32) 



where the coefficients of the interpolating function are given by 



+ A.J + B,ja,j) + ^h^a,jD^,j ] . (33) 
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